function  []=repeatability_study(input)
%  This code performs fits to RT dose response data in terms of tumor control
%  probability (TCP) using clinical information from a database: dose per 
%  fraction, number of RT fractions, treatment time, tumor control 
%  probability and number of patients/metastases treated. Two different TCP
%  formulations (Logit or Poisson) and several expressions of the surviving 
%  fraction (SF) models can be used with different functional dependencies (dep_a, dep_b):

%  log (SF)=n.*d.*alpha.*(1+a.*d.^dep_a+d.*(1+b.*d.^dep_b)/alphabeta),2)-lambda*max(0,t_trat-tk)

%  The fitting optimization process is n repeated times for each SF model. 

%  Input: One single string with input parameter values separated by spaces:
%  1) type of model (Logit or Poiss for logistic of Poisson TCP formulations,
%  respectively), 2) treatment  site (brain or lung_), 3) gamma or g value for
%  Logit or Poisson models, 4) alphabeta in Gy, 5) tumor accelerated repopulation 
%  kickoff time in days, 6) number of iterations in the simulated annealing  optimization,
%  8) simulated annealing cooling rate, and 9) n, number of desired repeated optimizations, 
%  multiple of 10.


%  Input examples: 'Poiss brain 3 10 28 60 0.965 20', 'Logit lung_ 0.83 10 28 100 0.975 20'


%  The data.mat file is loaded containing two clinical data arrays data_TCPexp_lung and data_TCPexp_brain
%  with a NSCLC and a brain metastases RT treatment cohorts dose response databases respectively.
%  Two additional vectors (dep_a and dep_b) specify the functional  dependencies of 26 SF models used.


%  Output: Matlab structure with:
 
%  1.	results_array: Results summary array with n x 26 bestfitting AICc values plus 1x 26 model's
%       AICc relative standard deviation, 1x 26 model's minimum AICc value from n repeated
%       simulations, 1x26 wi values associated to the best fits from n repeated simulations,
%       1x26 mean wi values calculated for each model, 1x26 wi standard deviation values for 
%       each model, 1x26 EVRs  associated to best fits from n repeated simulations, 1x26 mean
%       EVR values calculated for each model from n repeated simulations, and 1x26 EVR standard 
%       deviation for each model from n repeated simulations.

%   2.  n Model.sim_n Matlab structures with 8 fields each:

%       2.1	bp: 26 model's parameter values: 
%           - Poisson TCP model: bp=[alpha alphabeta a b g lambda kickoff_time TCPmax] 
%             Optimized parameters: alpha, lambda and a and b (when dep_a and dep_b are not zero, respectively)
%           - Logistic TCP model: bp=[alphabeta a b D50 gamma lambda kickoff_time TCPmax]
%             Optimized parameters: D50, lambda and a and b (when dep_a and dep_b are not zero, respectively)
%       2.2	Loglike: 26 model's best fitting loglikelihood values. 
%       2.3	aic: 26 model's best fitting AICc values
%       2.4	prob: fitting probability associated to each point of the database for the 26 models  
%       2.5	TCPmodel: Best fitting TCP model value associated to each point of the database (26 models)
%       2.6	EQD2: EQD2 associated to the each point of the database for the 26 models
%       2.7	Effect: Treatment effect, Log(SF),associated to the each point of the database for the 26 models.
%       2.8	Evol_loglike: 26 vectors of increasingly bestfitting loglike values calculated during the fitting optimization process.


tic;
%Database load
load('data.mat')
%Input parameters definition
inputs=strsplit(input);
model=inputs(1); %'Poiss' or 'Logit'
treatment_type=inputs(2); % 'brain' or 'lung_'
gamma=inputs(3);
alphabeta=inputs(4); %alpha/beta LQ model parameters ratio
tk=inputs(5); %Tumor repopulation kickoff time 
n_loop=str2double(inputs(6)); %simulated annealing iterations (100 for Logit, 60 for Poisson used in this work)
dT=str2double(inputs(7)); %simulated annealing cooling rate (0.975 for Logit, 0.965 for Poisson used in this work)
n=str2double(inputs(8)); %n� of optimizations (multiples of 10 used)
%Output file name
name=strcat('repeatability_', model, '_', treatment_type, '_g_', gamma, '_alphabeta_',alphabeta,'_tk_', tk, '.mat'); 
alphabeta=str2double(alphabeta);
gamma=str2double(gamma);
tk=str2double(tk);

if treatment_type{1}=='brain'
    data_TCPexp=data_TCPexp_brain;
    if model{1}== 'Logit' 
        %Initial parameter values for the TCP model: [alphabeta a b D50 gamma lambda kickoff_time TCPmax]
        param=[alphabeta -1e-5 1e-5 25 gamma log(2)/6 tk 0.95];
    elseif  model{1}=='Poiss'
        %Initial parameter values for the Poisson TCP model: [alpha alphabeta a b g lambda kickoff_time TCPmax]
        param=[0.36 alphabeta -1e-5 1e-5 gamma log(2)/6 tk 0.95]; 
    end
    
elseif treatment_type{1}=='lung_'
    data_TCPexp=data_TCPexp_lung;
    if model{1}=='Logit'
        %Initial parameter values for the TCP model: [alphabeta a b D50 gamma lambda kickoff_time TCPmax]
        param=[alphabeta -1e-5 1e-5 60 gamma log(2)/6 tk 0.95];
    elseif  model{1}=='Poiss'
        %Initial parameter values for the TCP model: [alpha alphabeta a b g lambda kickoff_time TCPmax]
        param=[0.18 alphabeta -1e-5 1e-5 gamma log(2)/6 tk 0.95];
    end
end

%ind_opt = indices of the param vector corresponding to parameters free to change during the optimization process.
    if model{1}=='Logit'
        logistic=1; 
%Optimized parameters: alpha, lambda and a and b (when dep_a and dep_b are not zero, respectively)
        ind_opt=[2 3 4 6];
    elseif model{1}=='Poiss'
        logistic=0;
%Optimized parameters: D50, lambda and a and b (when dep_a and dep_b are not zero, respectively)
        ind_opt=[1 3 4 6];
    end

%Initial simulated annealing temperature
T0=1;
     for i=1:n
        sim=['sim' num2str(i)];
        [Model.(sim).bp,Model.(sim).Loglike,Model.(sim).aic,Model.(sim).prob,Model.(sim).TCPmodel,Model.(sim).EQD2,Model.(sim).Effect,Model.(sim).evol_loglike]=MaxLikelihood_TCP(data_TCPexp,param,ind_opt, T0, dT, n_loop,logistic,dep_a,dep_b,treatment_type{1});
     end
%AIC array with nx26 AIC values corresponding to the bestfitting models
%calculated in each optimization
Model.aic=[Model.sim1.aic; Model.sim2.aic; Model.sim3.aic; Model.sim4.aic; Model.sim5.aic; Model.sim6.aic; Model.sim7.aic; Model.sim8.aic; Model.sim9.aic; Model.sim10.aic];
      if n/10>1
        for i=1:(n/10-1)
            extra=[Model.(['sim' num2str(i) '1']).aic; Model.(['sim' num2str(i) '2']).aic; Model.(['sim' num2str(i) '3']).aic; Model.(['sim' num2str(i) '4']).aic; Model.(['sim' num2str(i) '5']).aic; Model.(['sim' num2str(i) '6']).aic; Model.(['sim' num2str(i) '7']).aic; Model.(['sim' num2str(i) '8']).aic; Model.(['sim' num2str(i) '9']).aic; Model.(['sim' num2str(i+1) '0']).aic];
            Model.aic=[Model.aic; extra];
        end
      end
%Calculation of 26 Akaike weights from bestfitting models amongst n repeated optimizations
[Model.wimin]=weight_calculator(min(Model.aic));
%Calculation of nx26 Akaike weights from n optimizations
[Model.wi]=weight_calculator(Model.aic);
 
 for i=1:n; Model.ev(i,:)=Model.wi(i,:)./Model.wi(i,1); end
%Definition of results array:
 Model.results_array=[Model.aic; std(Model.aic)./mean(Model.aic); min(Model.aic); Model.wimin; mean(Model.wi); std(Model.wi); Model.wimin./Model.wimin(1); mean(Model.ev); std(Model.ev) ];
 Model=rmfield(Model, 'aic'); Model=rmfield(Model, 'wimin'); Model=rmfield(Model, 'ev'); Model=rmfield(Model, 'wi');
 time=toc;
 %Output file saved
 save(name{1},'Model','time')
end
